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We apply a hybrid Molecular Dynamics and mesoscopic simulation technique to study the steady- 
state sedimentation of hard sphere particles for Peclet numbers (Pe) ranging from 0.08 to 12. Hy- 
drodynamic back-flow causes a reduction of the average sedimentation velocity relative to the Stokes 
velocity. We find that this effect is independent of Pe number. Velocity fluctuations show the ex- 
pected effects of thermal fluctuations at short correlation times. At longer times, non-equilibrium 
hydrodynamic fluctuations are visible, and their character appears to be independent of the ther- 
mal fluctuations. The hydrodynamic fluctuations dominate the diffusive behavior even for modest 
Pe number, while conversely the short-time fluctuations are dominated by thermal effects for a 
surprisingly large Pe numbers. Inspired by recent experiments, we also study finite sedimentation 
in a horizontal planar slit. In our simulations distinct lateral patterns emerge, in agreement with 
observations in the experiments. 

PACS numbers: 05.40.-a,82.70.Dd,47.11.-j,47.20.Bp 



I. INTRODUCTION 

The steady-state sedimentation of spheres in a viscous 
medium at low Reynolds (Re) number is an important 
model problem in non-equilibrium statistical mechanics, 
exhibiting subtle and interesting physics [E El 9 ■ Some 
properties are relatively straightforward to determine. 
For example, the sedimentation velocity Vg of a single 
sphere was first calculated over 150 years ago by George 
Gabriel Stokes [1] to be V§ — fga 2 (p c — p) /v> where a 
and p c are the radius the density of the sphere, g is the 
gravitational acceleration, and p and rj are the density 
and viscosity of the fluid. On the other hand, even the 
first order effect of finite volume fraction <f> = |7rn c a 3 (n c 
is the particle number density), was not calculated until 
1972 when Batchelor [j| showed that 



V s = (1 - 6.550 + (C(j) 2 )) . 



(1) 



The effect of a finite volume fraction on sedimentation is 
dominated by long-ranged hydrodynamic forces that de- 
cay with interparticle distance r as slowly as r _1 . These 
forces are hard to treat analytically because they can eas- 
ily lead to spurious divergences. Eq. (Q]) also highlights 
the strong effect of the hydrodynamic forces. For exam- 
ple, a naive application of this lowest order result would 
suggest that all sedimentation should stop at 4> ~ 0.15. 
Of course this is not true since there are important higher 
order corrections in (j) whose calculation remains an ac- 
tive topic of research [6[- 

If the influence of hydrodynamics on the average sedi- 
mentation velocity at finite volume fraction is non-trivial 
to calculate, then the fluctuations around that average 
would appear even more formidable to determine. In 
a remarkable paper, Caflish and Luke Q used a simple 



scaling argument to predict that, for a homogeneous sus- 
pension, the velocity fluctuations 5V = V — Vg should 
diverge as (SV 2 ) ~ L, where L is the smallest container 
size. This surprising result stimulated much theoretical 
and experimental work, as well as no small amount of 
controversy Particle velocimetry experiments clearly 
show the existence of large-scale velocity fluctuations, 
which manifest as "swirls" [I S, E9, El EI- The ex ~ 
periments of Nicolai et al. [1| and Segre et al. @ suggest 
that while for small containers the velocity fluctuations 
do indeed grow linearly in L, for larger containers the 
velocity fluctuations saturate (see however [11 EH ) . The 
reasons (if any) that this should be observed have been 
the subject of sustained theoretical debate. It was shown 
by Koch and Shaqfeh [3] that hydrodynamic interac- 
tions can be screened if the colloids exhibit certain long- 
ranged correlations reminiscent of those found for electro- 
static systems. A number of theories have been proposed 
to generate such correlations in the bulk, including a cou- 
pled convective-diffusion model by Levine et al. [151 ] that 
generates a noise-induced phase-transition to a screened 
phase. Another class of theories focuses on the container 
walls. For example Hinch [l(| has argued that the bottom 
of a vessel will act as a sink for fluctuations, a prediction 
that appears to be confirmed by computer simulations 
[H EH • Other authors have emp hasized the importance 
of stratification [12, El, H3, HI HI and polydispersity 

mm. 

Most of the theoretical studies of sedimentation de- 
scribed above have focused on the non-Brownian limit 
where thermal fluctuations are negligible. This can be 
quantified by defining the Peclet number (Pe) 
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Via 



Mtga 
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where Dq is the equilibrium self-diffusion constant and 
Mb = ^ira 3 (p c — p) is the particle's buoyant mass. The 
non-Brownian limit then corresponds to Pe = oo. Be- 
cause Pe scales as (p c — p) a 4 , the very large Pe numbers 
needed to approximate the non-Brownian limit are easily 
achieved by increasing particle size. 

The Pe number is directly related to the gravitational 
length l g — kBT/(Mbg) = a/Ye. For this reason, the 
criterion Pe < 1 is often used to define the colloidal 
regime since, roughly speaking, one would expect from 
the barometric law that particles would then be dispersed 
throughout the solution. For example, for polystyrene 
spheres in water Pc = 0(1) for a « 1/im. In experi- 
ments, the density difference Ap = p c — p can be adjusted 
by density matching so that the Pe number can also be 
tuned quite accurately for a given a. 

In contrast to most previous theoretical and computa- 
tional studies, which have focused on the non-Brownian 
Pe = co limit, in this paper we study steady state sedi- 
mentation at the moderate Pe numbers relevant for the 
colloidal regime. In this regime the particles experience 
both (random) thermal fluctuations and (deterministic) 
hydrodynamic fluctuations, and a key question will be 
how these two kinds of fluctuations interact. 

We employ Stochastic Rotation Dynamics (SRD) 0, 
l25l |26| to describe the solvent, and a Molecular Dynam- 
ics scheme to propagate the colloids. Such a hybrid tech- 
nique was first employed by Malevanets and Kapral [2^ | , 
and recently used to study colloidal sedimentation by 
ourselves [28[ and by Hecht et al. [2{|. In section HT1 we 
briefly recap the salient details of our simulation method. 

In section IIIII we study the average sedimentation ve- 
locity. Our principle finding is that this follows exactly 
the same trend with volume fraction <fi as found for the 
Pe — 00 non-Brownian limit. In other words, the ef- 
fects of backflow are completely dominated by the hy- 
drodynamic interactions (HI), even when the Brownian 
forces are, on average, much stronger. In sections IIVI 
and [V] we investigate in some detail the velocity fluctua- 
tions (SV 2 ). We find that the thermal and hydrodynamic 
fluctuations appear to act independently of each other. 
Their effects are additive, at least in the accessed simu- 
lation regime, where the hydrodynamic fluctuations are 
unscreened. Some of these results have appeared earlier 
in a short letter [28j |. but here they are treated in much 
more detail. In section IVT1 we calculate the self-diffusion 
coefficient, highlighting the effects of hydrodynamic dis- 
persion. In section IVlTl we briefly consider the case of fi- 
nite sedimentation in a horizontal planar slit. We observe 
distinct lateral patterns, in agreement with recent laser 
scanning confocal microscopy. In section I Villi we dis- 
cuss the importance of thermal fluctuations over hydro- 
dynamic fluctuations. Finally, in section Hxl we present 
our conclusions. 



II. HYBRID MD-SRD COARSE-GRAINED 
SIMULATION METHOD 

The time- and length-scale differences between col- 
loidal and solvent particles are enormous: a typical col- 
loid of diameter 1 pm will displace on the order of 10 10 
water molecules. Clearly, some form of coarse- graining 
of the solvent is necessary. In this paper we use SRD 
to efficiently describe the dynamics of the solvent. The 
colloids are coupled to the solvent through explicit in- 
teraction potentials. We have recently performed an ex- 
tensive validation of this method 26] . We will therefore 
only reproduce the most important conclusions. 



A. Solvent-solvent interactions 

In SRD, the solvent is represented by a large number 
Nj of particles of mass mf. We will call these fluid parti- 
cles, with the caveat that, however tempting, they should 
not be viewed as some kind of composite particles or clus- 
ters made up of the underlying molecular fluid. The par- 
ticles are merely a convenient computational device to 
facilitate the coarse-graining of the fluid properties [26| . 

In the first step, the positions and velocities of the fluid 
particles are propagated by integrating Newton's equa- 
tions of motion. The forces on the fluid particles are 
generated by external forces generated by gravity, walls, 
or colloids. Direct forces between pairs of fluid parti- 
cles are, however, excluded; this is the main reason for 
the efficiency of the method. After propagating the fluid 
particles for a time At c , the second step of the algorithm 
simulates the collisions between fluid particles. The sys- 
tem is partitioned into cubic cells of volume ag. The 
velocities relative to the center of mass velocity v cm of 
each separate cell are then rotated: 

Vj H-> Vcm + R (v* - v cm ) . (3) 

R is a rotation matrix which rotates velocities by a fixed 
angle a around a randomly oriented axis. The aim of the 
collision step is to transfer momentum between the fluid 
particles. The rotation procedure can thus be viewed 
as a coarse-graining of particle collisions over time and 
space. Because mass, momentum, and energy are con- 
served locally, the correct (Navier-Stokes) hydrodynamic 
equations are captured in the continuum limit, including 
the effect of thermal noise [3] ■ 

Ihle and Kroll [IE] pointed out that at low temper- 
atures or small collision times At c the transport coef- 
ficients of SRD show anomalies. These anomalies are 
caused by the fact that fluid particles in a given cell can 
remain in that cell and participate in several collision 
steps. They showed that under these circumstances the 
assumption of molecular chaos and Galilean invariancc 
are incorrect. They also showed how the anomaly can 
be entirely cured by applying a random shift of the cell 
coordinates before the collision step. It is then possible 
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to analytically calculate the shear viscosity of the SRD 
fluid [30(. Such expressions are very useful because they 
enable us to efficiently tune the viscosity of the fluid, 
without the need of trial and error simulations. 



B. Colloid-colloid and colloid-solvent interactions 

In the simulation, colloidal spheres of mass M are 
propagated through the Velocity Verlet algorithm [3l[ 
with a time step A£md- The colloids are embedded in 
the fluid, and interact with the fluid particles through a 
repulsive (Weeks-Chandler- Andersen) potential: 



<Pcf(r) 




(r < 2 VV C/ ) 
(r > 2 VV C/ ) 



(4) 

The colloid-colloid interaction is represented by a similar, 
but steeper, repulsive potential: 




,48 



. 24 



(r < 2 1 /24 (7cc) 
(r > 2 1 /24 C7cc ) 



(5) 

The latter potential makes the colloids hard enough to 
approximate hard spheres, yet smooth enough to enable 
accurate integration of the equations of motion with a 
time step AtMD close to the collision time interval At c . 

Because the surface of a colloid is never perfectly 
smooth, collisions with fluid particles transfer angular 
as well as linear momentum. These interactions may be 
approximated by stick boundary conditions. We have 
studied several implementations of stick boundary con- 
ditions for spherical colloids [13] and derived a version of 
stochastic boundary conditions which reproduce linear 
and angular momentum correlation functions that agree 
with Enskog theory for short times and hydrodynamic 
mode-coupling theory for long times. Nevertheless, in 
this paper we use the radial interactions described in 
Eq. These do not transfer angular momentum to 

a spherical colloid and so induce effective slip boundary 
conditions. For many of the hydrodynamic effects we will 
discuss here the difference with stick boundary conditions 
is quantitative, not qualitative, and also well understood. 
For example, we have confirmed that the flowfield around 
a single sedimenting sphere decays, to first order, like 
a/(2r) for a slip boundary sphere [261 ] . whereas it decays 
like 3a/ (4r) for a stick boundary sphere. 

To avoid (uncontrolled) depletion forces, we routinely 
choose the colloid-fluid interaction range <r c / slightly be- 
low half the colloid diameter a cc /2 [26j . There is no a- 
priori reason why the hydrodynamic radius should be 
exactly half the particle-particle hard-core diameter for 
a physical system. For charged systems, for example, the 
difference may be substantial. An additional advantage 
of this choice is that more fluid particles will fit in the 
space between two colloids, and consequently lubrication 
forces will be more accurately represented. We have con- 



firmed that with our parameters SRD resolves the ana- 
lytically known lubrication forces down to gap widths as 
small as a/5. The agreement at small distances is caused 
also by repetitive collisions of the fluid particles trapped 
between the two surfaces. But at some point the lubrica- 
tion force will break down. An explicit correction could 
be applied to correctly resolve these forces for very small 
distances, as was implemented by Nguyen and Ladd [33| 
for Lattice Boltzmann dynamics. However, in this paper 
our choice of a c f is small enough for SRD to sufficiently 
resolve lubrication forces up to the point where the direct 
colloid-colloid interactions start to dominate 12611. 



C. Time scales and hydrodynamic numbers 

Many different time scales govern the physics of a col- 
loid of mass M embedded in a solvent. Hydrodynamic 
interactions propagate by momentum diffusion and also 
by sound. The sonic time is the time it takes a sound 
wave to travel the radius of a colloid, t cs — a/c s , where 
c s is the speed of sound. The kinematic time, on the 
other hand, is the time it takes momentum to diffuse 
over the radius of a colloid, t v = a 2 fv, where v is the 
kinematic viscosity of the solvent. For a colloid of radius 
a = 1/im in water, r cs 10~ 9 s and t u ~ 10~ 6 s. 

The next time scale is the Brownian time tb = M/£s, 
where £g = Qirrja is the Stokes friction for stick bound- 
ary conditions, or A-n-qa for slip boundary conditions. It 
measures the time for a colloid to lose memory of its ve- 
locity (see, however, (26|). The most relevant time scale 
for Brownian motion is the diffusion time td = a /Do, 
which measures how long it takes for a colloid to diffuse 
over a distance a in the absence of flow. For a colloid of 
a = Ifim in water, td ~ 5 s. 

When studying sedimentation, the Stokes time is the 
time it takes a single colloid to advect over its own radius, 
ts = a /Vs- The Stokes time and the diffusion time are 
related by the Peclet number: Pc = Tofts- If Pc 3> 1, 
then the colloid moves convectively over a distance much 
larger than its radius a in the time td that it diffused 
over the same distance. For Pe <C 1, on the other hand, 
the opposite is the case, and the main transport mech- 
anism is diffusive. It is sometimes thought that for low 
Pe numbers hydrodynamic effects can safely be ignored, 
but this is not always true, as we will show. 

In summary, in colloidal suspensions we encounter a 
range of time scales, ordered like t cs < tb < t„ < 
( T D,ts), where ts may be smaller or larger than td de- 
pending on Pe, and where we have assumed p c ~ p to 
justify tb < T„. The entire range of time scales can 
span more than 10 orders of magnitude. Thankfully, it is 
not necessary to exactly reproduce each of the different 
time scales in order to achieve a correct coarse- grai ning 
of colloidal dynamics. We can "telescope down" [26| the 
relevant time scales to a hierarchy which is compacted 
to maximize simulation efficiency, but sufficiently sepa- 
rated to correctly resolve the underlying physical behav- 
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ior. Keeping the relevant time-scales an separated by 
about an order of magnitude should suffice. 

Similar arguments can be made for various hydrody- 
namic numbers. For example, the Re number of sedi- 
menting colloidal particles is normally very low, on the 
order of 10 -5 or less. But there is no need to take such a 
low value since many relative deviations from the zero-Re 
Stokes regime scale with Re 2 . Exactly how big an error 
one makes depends on what one is investigating, but for 
our purposes we will take Re < 0.4 as an upper bound. 
We have shown (2(| that for the friction on a sphere in- 
ertial effects are unimportant up to Re ss 1. We note 
that our upper bound on Re also ensures that the time 
hierarchy condition t„ < (is, 175) is fulfilled. In principle 
Pe can be whatever we like as long as Re remains low 
and the hierarchy is obeyed. 

In our simulations we choose an average number of 
fluid particles per collision volume equal to 7 = 5, a col- 
lision interval At c = 0.1 (in units of to = ao^/mj /fcgT), 
and a rotation angle a = tt/2, leading to a kinematic 
viscosity v = 0.5 a^/t - We choose a colloidal mass 
M = 125mj, and interaction parameters a c f — 2ao, 
a cc — 4.3a , and e = 2.bk B T. We have verified that this 
choice leads to a small relative error in the full velocity 
field, and that we can quantitatively calculate the the 
observed friction on a colloid [26] . Note that this friction 
is somewhat lower than expected on the basis of a hydro- 
dynamic radius set equal to er c / = 2a,Q. This is due to ad- 
ditional Enskog friction effects, where the different con- 
tributions to the friction add "in parallel" , as explained 
in Ref . [26[ . The resulting effective hydrodynamic radius 
a = 1.55ao will be used throughout this paper. The time 
scales in our simulations are well separated: t cs = 1.2io, 
t b = 2.5to, Ty = 4.8*0, and t d = 120i . 



III. AVERAGE SEDIMENTATION VELOCITY 

Sedimentation simulations were performed in a peri- 
odic box of dimensions L x = L y = 32ao and L z = 96ao 
(approx. 21 x 21 x 62a), with N = 8 to 800 colloids. 
The number of SRD particles was adjusted so that the 
free volume outside the colloids contained an average of 
5 particles per coarse-graining cell volume a§. This cor- 
responds to a maximum of Nf w5x 10 5 SRD particles. 
A gravitational field g, applied to the colloids in the z di- 
rection, was varied to produce different Peclet numbers, 
ranging from Pe = 0.08 to Pe = 12. At the same time 
the Reynolds numbers ranged from Re = 0.003 to 0.4. 
Note that the absence of a wall at the top and bottom of 
the simulation box necessitates an additional constraint 
to keep the center-of-mass of the entire system fixed. 

Right after the simulations start, the colloidal positions 
and velocities have not yet acquired their steady-state 
distributions. We monitored block averages (in time) of 
the sedimentation velocity and the behavior of sedimen- 
tation velocity fluctuations, which will be discussed in 
the next section. We verified that there was no drift in 



10 rv---'- 
0.8 - 



0.6 - 



> 

> m 0.4 



0.2 - 



0.0 



X 


Pe 


= 0.08, L/a = 20.6 


* 


Pe 


= 0.2 


V 


Pe 


= 0.8 


O 


Pe 


= 4 


O 


Pe 


=12 


A 


Pe 


=12, L/a = 31.0 


□ 


Pe 


=12, L/a = 41 .3 



0.00 



0.05 



0.10 



0.15 



FIG. 1: Average sedimentation velocity, Vs normalized by 
the Stokes velocity V$, as a function of volume fraction (f> 
for various Peclet numbers and system sizes. Dashed lines 
correspond to the semi-empirical Richardson- Zaki law (1 — 
(j>) n , with n = 4.7 for the upper and n = 6.55 for the lower 
line. The dotted line is another theoretical prediction taking 
higher order HI into account [6|. Ignoring hydrodynamics 
leads to Vs/Vg =1 — (dash-dotted line). 



these properties after about 100 Stokes times ts, corre- 
sponding to sedimentation down the height of about two 
periodic boxes. The absence of any drift indicated that 
the suspensions were now in steady state. 

The simulations were subsequently run between 200 ts 
for Pe = 0.08 to 30,000 t s for Pe = 12. To check that 
our system is large enough, we performed some runs for 
1.5 and 2 times the box size described above, finding no 
significant changes in our conclusions. 

The average sedimentation velocity Vs for different 
Peclet numbers and system sizes as a function of hydro- 
dynamic packing fraction cf> = ^nn c a 3 is shown in Fig.Q] 
The results are normalized by the Stokes velocity Vg (the 
sedimentation velocity of a single particle in the simu- 
lation box), resulting in the so-called hindered settling 
function. At low densities the results are consistent with 
the result found by Batchelor while at higher den- 
sities they compare well with a number of other forms 
derived for the Pe — > 00 limit. In most experiments the 
hindered settling function is well described by the semi- 
empirical Richardson-Zaki law VsJVS = (1 — <fi) n , with 
n ranging between 4.7 and 6.55 [l|, y|. Our results fall 
between these two extremes. The results compare par- 
ticularly well with a theoretical prediction by Hayakawa 
and Ichiki [(|, taking higher order hydrodynamic inter- 
actions into account. 

One might naively expect that the effect of HI becomes 
weaker for Pe < 1. Taking into account only Brownian 
forces would result in Vs = Vg(l — (f>) (because of flux 
conservation) , which heavily underestimates backflow ef- 
fects. However, we observe that the results for all Peclet 
numbers 0.08 < Pe < 12 lie on the same curve. We em- 
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FIG. 2: Main plot: colloid radial distribution function g(r) for 
4> = 0.04 at low (0.2) and high (12) Peclet number. There is no 
significant difference between the two g(r)'s. Inset: structure 
factor for the same systems. At Pe = 12, small deviations 
are found for perpendicular (open circles) and parallel (closed 
circles) wave vectors. 



phasise that these results are normalised by the Stokes 
velocity Vg of a single sphere, which itself decreases with 
decreasing Peclet number. The important point is that 
the additional hindrance caused by hydrodynamic inter- 
actions is observed to be unaffected by the actual Pe 
number. A reason for this could be that the average sed- 
imentation velocity is determined predominantly by the 
(time-averaged) distribution of distances between the col- 
loids. If this is so, then the particle motion generated by 
the external field must not lead to a significant change 
in the microstructure. That this is indeed the case is 
shown in Fig. [2 where the main plot shows the colloidal 
radial distribution function at volume fraction <f) = 0.04 
for Peclet numbers 0.2 (stars) and 12 (circles). For Pe 
= 0.2 the result is indistinguishable from equilibrium re- 
sults, and for Pe = 12, despite the fact that the external 
field is quite strong, the average number of neighbour- 
ing particles at a certain distance from a specific particle 
changes only very slightly as compared to equilibrium. 
The inset of Fig. [2] shows the structure factor for the 
same system. At Pe = 12, small deviations are found for 
perpendicular (open circles) and parallel (closed circles) 
wave vectors, but again the differences are not very large. 
Here we already note that all of these systems are in the 
unscreened regime. 



IV. SPATIAL CORRELATIONS IN 
FLUCTUATIONS 

We next discuss velocity fluctuations around the av- 
erage. In colloidal systems the instantaneous velocity 
fluctuations <5V = V — V s are dominated by thermal 
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FIG. 3: Spatial correlation functions of the parallel (z) com- 
ponent of the velocity fluctuations as a function of distance 
perpendicular (a) and parallel (b) to the external field, for 
three different volume fractions (0 = 0.02 (grey symbols), 
cj> = 0.04 (white), (j> = 0.086 (black)) and different Peclet 
numbers. The correlation functions are scaled with Vg to 
emphasize hydrodynamic fluctuations. The insets show how 
C z (r), scaled with fcsT/M, increases with Pe. 



fluctuations, with a magnitude determined by equiparti- 
tion: 



AV£ = k B T/M. 



(6) 



To disentangle the hydrodynamic fluctuations from ther- 
mal fluctuations, we describe spatial and temporal cor- 
relations in the velocity fluctuations. The spatial corre- 
lation of the z component (parallel to the sedimentation 
direction) of the velocity fluctuations can be defined as 



C z (r) ee (SV z (0)SV z (r)) 



(7) 



where (...) represents an average over time and over all 
colloids. The distance vector r is taken perpendicular 
to sedimentation, C z (x), or parallel to it, C z (z). Note 
that we will not normalize the correlation functions by 
their initial values. Rather, we will normalize them by 
values which have a more physical meaning, such as the 
squared sedimentation velocity Vg, or the thermal fluc- 
tuation strength ksT/M. 

In Fig. [3] we plot C z (r) , which shows a positive spa- 
tial correlation along the direction of flow, and an anti- 
correlation perpendicular to the flow, very much like that 
observed in the experiments of Nicolai et al. [|[ . The in- 
set of Fig. E3a) shows that at Pe = 0.8 the correlation in 
the perpendicular direction, C z (x), is almost negligible 
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FIG. 4: Spatial correlation functions of the parallel (z) com- 
ponent of the velocity fluctuations as a function of distance 
perpendicular to the external field, for different system sizes. 
In the inset, distance is scaled with the horizontal box size L. 
All simulations were performed at <j> = 0.04 and Pe = 8. 



compared with the thermal fluctuation strength fc^T/M, 
whereas for larger Pe, distinct regions of negative ampli- 
tude emerge, which grow with increasing Pe. Similarly, 
the inset of Fig. [3](b) shows correlations in the parallel 
direction that rapidly increase with Pe. For the highest 
Peclet numbers studied (4 < Pe < 12), the amplitudes of 
these correlations grow proportionally to Vf , as shown in 
the main plots of Fig. [3J Unfortunately, because the divi- 
sion by Vg amplifies the statistical noise, we are unable to 
verify whether this scaling persists for Pe < 4. The min- 
imum in Fig. 02a) is at about half the box width (this is 
also the reason why no data points could be collected for 
x > 15a). This suggests that the velocity fluctuations are 
unscreened and only limited by our box dimensions (see 
[13|). We confirm this in Fig. where it is seen that the 
correlation length scales linearly with box dimensions. 



V. TEMPORAL CORRELATIONS IN 
FLUCTUATIONS 

Similarly to the spatial correlations of the previous sec- 
tion, the temporal correlation of the z component of the 
velocity fluctuations can be defined as 



C z {t) = (SV z (0)SV z (t)) , 



(8) 



where now t is a correlation time and (...) denotes an 
average over all colloids and all time origins. Fig.[5]shows 
the temporal correlation functions along the direction of 
sedimentation on a linear scale. Clearly the correlation 
is increasing with increasing Pe number. To investigate 
this in more detail, we plot the temporal correlation on 
a log-log and log-linear plot in Fig. [5] 
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FIG. 5: Temporal correlation functions of the z component 
of the velocity fluctuations for <f> = 0.02 and different Peclet 
numbers. Time is scaled with the Brownian relaxation time 
tb = M/£ and the velocities are scaled with the thermal 
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FIG. 6: Temporal correlation functions of the z component 
of the velocity fluctuations for <f> = 0.02 and different Peclet 
numbers, (a) Time is scaled with the Brownian relaxation 
time tb = M/£ and the velocities are scaled with the thermal 
fluctuation strength ksT /M. The straight line is the hydro- 
dynamic long time tail Bt~ s/2 with B" 1 = l2pk B T(-n-is)' i/2 
[371 ]. The results for Pe < 1 are indistinguishable, (b) Time is 
scaled with the Stokes time t$ = a/Vs and the velocities are 
scaled with Vg to highlight hydrodynamic velocity fluctua- 
tions. The straight line is a fit demonstrating the exponential 
decay of non-equilibrium hydrodynamic fluctuations. 



7 



At very short times the velocity de-correlation is quan- 
titatively described by Enskog dense-gas kinetic the- 
ory [111, [36| , which predicts the following decay: 



lhnC z (t) = AFf exp( 



where the Enskog friction coefficient is given by 



_ 8 {2nk B TMmf 
^ ~ 3 V M + m f 



1/2 



(9) 



(10) 



Eq. @ describes the velocity relaxation due to random 
collisions with the solvent particles. 

At intermediate times the temporal correlation follows 
the well known algebraic long time tail 



Clong(t) = Bt 



- R+-3/ 2 



(11) 



associated with the fact that momentum fluctuations dif- 
fuse away at a finite rate determined by the kinematic 
viscosity v. Analytical mode-coupling calculations yield 
aprefactorB- 1 = 12 pk B T '(tt jvf' 12 [33]. This exactly fits 
the low Pe (< 1) results in Fig. Eta) with no adjustable 
parameters. We note that similar agreement was found 
for the long-time tails for other parameter choices [HJ at 
equilibrium. Of course these simulations are all at finite 
Pe number, and so are out of equilibrium, but for small 
Pe the long-time tail dominates within the simulation 
accuracy that we obtain. 

In an experimental study on the sedimentation of non- 
Brownian (Pe — > oo) particles, Nicolai et al. @ found an 
exponential temporal relaxation of the form 

C*(t) = Ay£exp(-t/Tff). (12) 

This non-equilibrium hydro dynamic effect takes place 
over much longer time-scales than the initial exponential 
relaxation due to random collisions with the solvent par- 
ticles, i.e. th 3> M/£e- The double-logarithmic Fig. Eta) 
shows that a new mode of fluctuations becomes distin- 
guishable in our simulations for Pe > 1. In the log- linear 
Fig. Etb) the correlation functions are scaled with Vg to 
highlight the nonequilibrium hydrodynamic fluctuations. 
For Pe > 8 the fluctuations scale onto a single exponen- 
tial master curve, similar to the high-Pe experiments of 
Nicolai et al. 8|, whereas for lower Pe deviations are 
seen. From the exponential fit to Eq. l|12j) . we can es- 
timate the relaxation time tjj and the amplitude AV^ 
of the hydrodynamic fluctuations. These are shown in 
Fig. [7] for different volume fractions <j>, and in Fig. [8] for 
different box sizes L/a. The results are consistent with 
a scaling AVh/Vs oc y/L<j>/a and Tufts oc yJLj {<j>a). 

These scalings can be understood by a simple heuris- 
tic argument by Hinch et al. [381 ] akin to that used by 
Caflish and Luke 0: Suppose we consider the box vol- 
ume to consist of two equally large parts, each with a 
typical linear dimension of L. The average number of 
colloids in a volume of size L 3 is (N) — L 3 <f)/ (|7ra 3 ). Of 
course the colloids are free to move from one part to the 
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other; the division is entirely artificial. At low enough 
volume fraction <fi we assume that the colloidal positions 
are described by random Poisson statistics. The typical 
fluctuation in the number of particles will then be of or- 
der yj (N). The extra colloidal weight of order \J (N)Mbg 
in one part of the box causes this part to sediment faster 
than average. This is the hydrodynamic fluctuation re- 
ferred to before. The extra colloidal weight is balanced 
by the extra Stokes drag caused by the larger sedimenta- 
tion velocity, which is of the order of GtttiLAVh ■ Making 
use of Vs = Mbg / (Sirr/a) , we predict for the amplitude of 



the hydrodynamic fluctuations: 



L<j> 



(13) 



This is consistent with our observation. Of course the 
hydrodynamic fluctuation does not persist indefinitely. 
It will decorrelate on the order of the time needed to 
fall over its own length, for it will then encounter and 
mix with a region of average number density. The relax- 
ation time of the hydrodynamic fluctuation is therefore 
predicted to be 



L 2 



•h 



AV* 



^naL 2 

W4 



4-1 6 



a4> 



(14) 



where we have used ts = a /Vs. 

The above scaling argument does not fix the prefactors. 
Fitting with the data in Figs. [7] and [5] we find AVh ~ 
0.29V S [0(L/a)] 1/2 and t h w 0.33f s [cj>(a/ L)]~ 1/2 . It 
should be noted that the above results concern the ve- 
locity fluctuations parallel to the gravitational field (z). 
In a similar way we have estimated the perpendicular 
velocity fluctuations to be characterized by AVff iPorp « 
Q.lW s [^{L/a)] 1/2 and r^ porp « 0.15t s [<l>{a/ L)]' 1 ' 2 . 
Note that the ratio of parallel to perpendicular velocity 
fluctuations is approximately 1.8. This is in agreement 
with the experimental low (f> results on non-Brownian 
spheres by Nicolai et al. [|[ and by Segre et al. [§] , both of 
whom observed vertical fluctuations approximately twice 
the horizontal fluctuations in the same range of volume 
fractions. 



VI. DIFFUSION AND DISPERSION 

The equilibrium self-diffusion of a colloidal particle is 
related to its velocity correlation function through the 
following Green-Kubo equation: 

D (t)= [ (V x (r)V x (0))dT, (15) 
Jo 

where V x is a Cartesian component of the colloidal veloc- 
ity. For large enough times t the integral Do (t) converges 
to the equilibrium self-diffusion coefficient Do. 

During sedimentation, the diffusion is enhanced by the 
hydrodynamic fluctuations. In fact, the diffusion is no 
longer isotropic but tensorial. Focusing first on the com- 
ponent parallel to gravity, we define the parallel diffusion 
coefficient similarly to Eq. (fT5")) as the large time limit of 



D z (t) = / (SV z {T)8V z (0))dn 



(16) 



In Fig. [5]we show D z (t) normalized by the equilibrium 
value Do for a range of Pe numbers. Note that even 
though the hydrodynamic fluctuations may be small com- 
pared to C(0), they nevertheless have a significant con- 
tribution to the diffusivity because the time-scale th is 
much longer than t v . 
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FIG. 9: Time dependent self-diffusion coefficient parallel to 
gravity for different Pe numbers. Do is the equilibrium self- 
diffusion coefficient. 



To understand the total diffusivity, we make the fol- 
lowing addition approximation: 



D z = D n 



D 



H 



(17) 



where Do is equilibrium diffusion coefficient and Dh the 
dispersion due to non-equilibrium hydrodynamic fluctu- 
ations. The former can be approximated as a sum of 
Stokes and Enskog diffusion coefficients, see [26]. The 
non-equilibrium hydrodynamic dispersion can be esti- 
mated using the previous scaling arguments: 



D H w AV h t h oc Vsot/)^ 



L 



3/2 



(18) 



Taking the prefactors found in the previous section, and 
rewriting Vsa as PeDo, we therefore predict 



D z =D n 



1 + 0.03 Pe 



a/2 (L 



3/2' 



(19) 



for small enough, i.e. unscreened, systems. For small Pe 
(< 1) the self-diffusion coefficient is largely independent 
of Pe and equal to D , whereas for very large Pe (3> 1) 
it becomes proportional to Pe. This is confirmed in Fig. 
[TOl where the dashed lines show the Pe and <f> dependence 
of Eq. [H 

The diffusion in the plane perpendicular to gravity is 
also enhanced by the hydrodynamic fluctuations, similar 
to Eq. (p~9|) . but with a smaller prefactor of 0.004 instead 
of 0.03 (not shown). The ratio of hydrodynamic diffu- 
sivities, D/f/D^f iPcrp ps 7 is similar to what is found in 
the experiments of Nicolai et al. (8J for non-Brownian 
spheres. 

Although our simulations are in the unscreened regime, 
it is interesting to also consider the hydrodynamic contri- 
bution to the diffusion coefficient in the screened regime. 
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Pe 

FIG. 10: Self-diffusion coefficient parallel to gravity versus 
Peclet number for different concentrations (j>. Do is the equi- 
librium self-diffusion coefficient. Dashed lines are predictions 
from Eq. 1(15)1. 

If we apply the experimental fits of Segre et al. Q for 
AVff and the correlation length £, then the simple scal- 
ing arguments above suggest that 

D H (xPeD , (20) 

which is independent of <f>. The exact pre-factor is hard 
to determine in the screened regime. Nevertheless, an 
estimate can be made if we assume that th has the same 
pre-factor in the experiments as we find in our simula- 
tions. For example, if we replace L/2, which measures the 
location of the minimum of the perpendicular correlation 
functions, with £ pe rp; its value for the screened regime [9J, 
then we find D H>peip /D w 1.1 Pe. For Dff, P arallel/-Do we 
expect a pre-factor several times larger. In the screened 
regime the hydrodynamic contributions to the diffusion 
should dominate for Pe> 1. In practice however, we ex- 
pect that for many colloidal dispersions effects such as 
polydispersity[l8| may temper the size of the swirls, and 
thus reduce the hydrodynamic contribution to diffusion. 



VII. FINITE SEDIMENTATION IN A 
HORIZONTAL PLANAR SLIT 

Up to this point we have focused on steady-state sedi- 
mentation by applying periodic boundary conditions and 
giving the system enough time to overcome transient flow 
effects. 

One may wonder what happens if the particles are con- 
fined and are not given enough time to reach steady-state. 
Very recently, Royall et al. [3!| studied nonequilibrium 
sedimentation of colloids in a horizontal planar slit, at 
a Peclet number of order 1, using laser scanning confo- 
cal microscopy. Among other things, they measured the 



time evolution of the one-dimensional colloid density pro- 
file p(z,t), where the z-axis is normal to the horizontal 
plane. Two cases were considered. In the first case an ini- 
tially homogenised sample was allowed to sediment to the 
bottom of the capillary. Good agreement was found with 
a dynamical density functional theory (DDFT) calcula- 
tion that included a density-dependent mobility function. 
In the second case they considered an equilibrated sam- 
ple turned upside down so that the previous sediment 
suddenly finds itself at the top of the capillary. In this 
case sedimentation proceeds in an entirely different fash- 
ion. A strong finger-like inhomogeneity was observed, 
accompanied by maze-like lateral pattern formation. 

Inspired by these experiments, we set up a box of size 
180 x 180 x 60ao (116 x 116 x 39a), with periodic bound- 
aries in the x and y direction, and with walls at the top 
and bottom in the z direction. (This corresponds to a 
height of about 32a, close to the experimental value of 
36a.) We add 6500 colloids (4> ~ 0.06) and apply an 
external field upwards such that Pe = 4. After reaching 
the equilibrium distribution, at t — we suddenly reverse 
the field, again at Pe = 4. We observe a maze-like lateral 
pattern, Fig. 1111 which shows striking similarities to the 
experimental observations [39( . The characteristic length 
of the maze-like lateral pattern is approximately equal 
to the height of the slit. It has been suggested [3!| that 
there may be a relation between this phenomenon and the 
swirls observed in steady-state sedimentation, but also 
that the swirls are reminiscent of a Rayleigh- Taylor type 
instability in two layered liquids, with the steep initial 
density gradient resembling a (very diffuse) liquid-liquid 
interface. With the current data, we cannot conclusively 
determine the origin of this instability. Nevertheless it is 
gratifying that our simulations produce such similar, and 
nontrivial, behaviour as the experiments under similar 
conditions. This can be viewed as an additional valida- 
tion of our simulation model. 

In Fig. [TJ] we analyse the time evolution of the one- 
dimensional density profile p(z, t). The crystal-like layers 
at the top plate for t < 0, disappear and then reappear 
again at the bottom of the plate. It would be interest- 
ing to compare these results to calculations using DDFT. 
Since the latter technique does not explicitly contain any 
long-ranged hydrodynamics, one would expect it to have 
difficulty in reproducing the swirls observed in the sim- 
ulations and experiments. Nevertheless, because both 
the initial and final states are constrained by equilibrium 
statistical mechanics (for which DFT is very accurate), 
the one-body density p(z, t) may not be a very sensitive 
measure of the more complex dynamics that arise from 
hydrodynamics. 



VIII. DISCUSSION 

As seen in Fig. [|5J the short time velocity fluctuations 
are dominated by thermal fluctuations at all Peclet num- 
bers studied. The relative strength of the t = thermal 



10 




10 
8 

— 6 

-t— ' 

N 

^ 4 
2 




~" 1 

-- t = 
.-. t = 4x D 

-t = 8x D 

— t= 12x n 




ii 
ii 
ii 
ii 
ii 
ii 
ii 
ii 
ii 
ii 
n 
1 1 
1 1 
1 1 
1 1 



40 



FIG. 12: Time evolution (right to left) of the one-dimensional 
density profile for sedimentation in a horizontal slit, p is nor- 
malised such that it equals 1 for a homogeneously filled slit. 
The final state (on the left) closely resembles the initial state 
(on the right), but is not shown for clarity. 



FIG. 11: An equilibrated sediment in a planar slit is turned 
upside down and allowed to sediment at Pe = 4. Shown here 
are the horizontal sy-plane and the corresponding vertical 
xz-plane at 6 different times after field reversal. The dashed 
line indicates the height z where the snapshots of the a;y-plane 
are taken. A strong fingerlike inhomogeneity develops quickly, 
accompanied by maze-like lateral pattern formation. 



and hydrodynamic velocity fluctuations follows from sim- 
ple scaling relations. Using 



RePe 



(Vsa) 2 



(21) 



which follows from the definitions of Pe and Re, together 
with Eqs. © and (jTHJ), the following relationship between 
hydrodynamic and thermal fluctuations emerges: 



—Mr w a Re Pe < 



^Lpc 
a p 



(unscreened) , 



(22) 



where the simplifying assumption that M c w ^irp c a 3 , 
with a the hydrodynamic rather than the physical radius, 
was also made. The numerical pre-factor a is small and 
can be extracted from Fig. [8] to be a « 0.05 for fluctua- 
tions parallel to the flow, and a « 0.015 for fluctuations 
perpendicular to the flow. 

The above scaling holds for the unscreened regime; in 
the screened regime the ratio of Vh to Vr will be smaller. 
Consider, for example, the experimental results of Segre 
et al. Q . If we take their fits to the scaling of the parallel 
fluctuations in the screened regime, together with the 
estimates Re = 5 x 10~ 5 and Pe « 2000, the scaling 
becomes: 



AVI 



6 2/3 RePe (screened) 



(23) 



for flows in the parallel direction. This suggests that 
this ratio is small in the experiments, from 2 x 10~ 4 
for <f) = 10~ 4 to 0.02 for <p = 0.1. So despite the fact 
that the Pe number in these experiments appears to be 
high, there is no need for an effective gravitational "tem- 
perature" 10] to thcrmalisc: at short correlation times 
the usual thermal fluctuations are still dominant. How- 
ever, because the product RePe scales with quite a high 
power of a, as fast as a 7 , the ratio AV^ /AV^ will increase 
rapidly for larger particles and gravitational temperature 
will become essential for thermalisation. 

When comparing parallel and perpendicular compo- 
nents it is important to mention that in numerical works 
where thermal fluctuations are neglected very strong 
anisotropics in velocity fluctuations, hydrodynamic re- 
laxation times, and diffusivities are often found. For 
example Ladd finds Dh pC r P ~ 25 in his Lat- 
tice Boltzmann simulations. This was attributed to pe- 
riodic boundary conditions. However, we also use peri- 
odic boundary conditions and find results much closer 
to experimental results (a diffusivity ratio of ~ 7). We 
therefore conclude that thermal fluctuations reduce the 
anisotropy. This could be tested in Lattice Boltzmann 
simulations by adding fluctuating stress (4lT |. 



IX. CONCLUSION 

In conclusion, we have studied the interplay of hydro- 
dynamic and thermal fluctuations using a novel simula- 
tion technique. The two types of fluctuations appear to 
act independently, at least in the unscreened regime. We 
find that hydrodynamic interactions are important for 
the average sedimentation velocity for Peclet numbers as 
low as 0.08, whereas thermal fluctuations may remain im- 
portant up to very large Peclet numbers. Neither may be 
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ignored for a significant range of Peclet numbers. We also 
calculate the hydrodynamic contributions to the diffu- 
sion coefficient, and find that with increasing Pe number 
they rapidly become much larger than the equilibrium 
diffusion coefficient. As an additional test of the method 
we studied finite sedimentation in a horizontal slit, and 
found characteristic lateral patterns in agreement with 
recent experiments. 
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